
;;================================================================
;; Simple routine to make comparison plots for the output of the
;; mass-function routine.  This removes junk plotting code from the
;; production routine.


restore,'./masses/iter_mfn.sav',/ver ; CHANGE TO LOCAL AT SOME POINT!!!
nreal = n_elements(iter)

;; Run some basic stats
alpha   = mean(iter.alpha)
e_alpha = stddev(iter.alpha)
mmin    = mean(iter.mmin)
e_mmin  = stddev(iter.mmin)
mtot    = mean(iter.mtot)
e_mtot  = stddev(iter.mtot)

message,'Alpha: '+string(alpha,e_alpha,format="(F0.2,' +- ',F0.2)"),/inf
print,m4_stat(iter.alpha)
message,'M_min: '+string(mmin,e_mmin,format="(E0.2,' +- ',E0.2)"),/inf
print,m4_stat(iter.mmin)
message,'M_tot: '+string(mtot,e_mtot,format="(E0.2,' +- ',E0.2)"),/inf
print,m4_stat(iter.mtot)

myps,'./masses/plots/iter_compare.eps',xsize=13.5
multiplot_xm,[3,1],mpcharsize=1.0
yr = [1.6,3.0]
sunsym = cgSymbol('sun')


;;====================================================================
;; Panel 1
xr=[100,1200]
cgPlot,alog10(iter.mmin),iter.alpha,psym='filledcircle',xst=5,$
       xtit='M!dmin!n  [M'+sunsym+']',yr=yr,/yst,symsize=0.7,/nodata,$
       ytit='Power-Law Slope   '+cgSymbol('alpha'),charsize=1.0,$
       xr=alog10(xr)
;; cgOplot,iter.mmin,iter.alpha,psym='filledcircle',symsize=0.7
density = HIST_2D(alog10(iter.mmin),iter.alpha, $
                  min1=alog10(xr[0]),max1=alog10(xr[1]),$
                  min2=yr[0],max2=yr[1],bin1=0.01,bin2=0.02)
print,m4_stat(density)
cgLoadct,13
plotimage,density,range=set_plot_range(density),/over,xst=4,yst=4

density = smooth(density,5)


;; Figure the contours for 68.27%, 95.45%, and 99.73% confidence regions
level = derive_contour_levels(density,nreal)
print,'LEVEL: ',level,max(density)

cgContour,density,levels=level,thick=3,/over,c_color='bisque',$
          c_annot=['99.7%','95.5%','68.3%'],c_charsize=0.8




al_legend,/top,/left,/clear,$
          ['N!diter!n = '+string(n_elements(iter),format="(E0.1)")]
cgAxis,xaxis=0,xr=xr,/xlog,charsize=1.0,/xst,$
       /save,xtit='M!dmin!n  [M'+sunsym+']'
cgAxis,xaxis=0,/xst,xtickformat='blank_axis',color='wt1'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',color='wt1'
cgAxis,yaxis=0,yr=yr,/yst,ytickformat='blank_axis',/save,color='wt1'
cgAxis,yaxis=1,ytickformat='blank_axis',color='wt1',/yst

vline,/h,2.35,thick=3,color='crimson',/xlog

cgPlots,3.21d2,2.24,psym='filleduptriangle',color='pur2',symsize=1.3
cgPlots,3.21d2,2.24,psym='openuptriangle',color='black',symsize=1.3

multiplot

;;====================================================================
;; Panel 2
xr=[1.25,1.95]
cgPlot,iter.mtot/1.d5,iter.alpha,psym='filledcircle',$
       xtit='M!dtot!n  [x10!u5!n M'+sunsym+']',xminor=2,$
       charsize=1.0,/nodata,yr=yr,/yst,symsize=0.7,xr=xr,/xst
;; cgOplot,iter.mtot/1.d5,iter.alpha,psym='filledcircle',symsize=0.7
density = HIST_2D(iter.mtot/1.d5,iter.alpha, $
                  min1=xr[0],max1=xr[1],$
                  min2=yr[0],max2=yr[1],bin1=0.005,bin2=0.02)
print,m4_stat(density)
cgLoadct,13
plotimage,density,range=set_plot_range(density),/over,xst=4,yst=4

density = smooth(density,5)


;; Figure the contours for 68.27%, 95.45%, and 99.73% confidence regions
level = derive_contour_levels(density,nreal)
print,'LEVEL: ',level,max(density)

cgContour,density,levels=level,thick=3,/over,c_color='bisque',$
          c_annot=['99.7%','95.5%','68.3%'],c_charsize=0.8



cgAxis,xaxis=0,/xst,xtickformat='blank_axis',color='wt1',$
       xr=xr,xminor=2,/save
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',color='wt1',$
       xr=xr,xminor=2
cgAxis,yaxis=0,yr=yr,/yst,ytickformat='blank_axis',/save,color='wt1'
cgAxis,yaxis=1,ytickformat='blank_axis',color='wt1',/yst

vline,/h,2.35,thick=3,color='crimson'

cgPlots,1.59,2.24,psym='filleduptriangle',color='pur2',symsize=1.3
cgPlots,1.59,2.24,psym='openuptriangle',color='black',symsize=1.3

multiplot

;;====================================================================
;; Panel 3
xr = [375,435]
cgPlot,iter.ndist,iter.alpha,psym='filledcircle',xtit='N!ddist!n',$
       charsize=1.0,yr=yr,/yst,symsize=0.7,/nodata,xminor=5,$
       xr=xr,/xst
;; cgOplot,iter.ndist,iter.alpha,psym='filledcircle',symsize=0.7
density = HIST_2D(iter.ndist,iter.alpha, $
                  min1=xr[0],max1=xr[1],$
                  min2=yr[0],max2=yr[1],bin1=1,bin2=0.02)
print,m4_stat(density)
cgLoadct,13
plotimage,density,range=set_plot_range(density),/over,xst=4,yst=4

density = smooth(density,5)


;; Figure the contours for 68.27%, 95.45%, and 99.73% confidence regions
level = derive_contour_levels(density,nreal)
print,'LEVEL: ',level,max(density)

cgContour,density,levels=level,thick=3,/over,c_color='bisque',$
          c_annot=['99.7%','95.5%','68.3%'],c_charsize=0.8



cgAxis,xaxis=0,/xst,xtickformat='blank_axis',color='wt1',$
       xr=xr,xminor=5,/save
cgAxis,xaxis=1,/xst,xtickformat='blank_axis',color='wt1',$
       xr=xr,xminor=5
cgAxis,yaxis=0,yr=yr,/yst,ytickformat='blank_axis',/save,color='wt1'
cgAxis,yaxis=1,ytickformat='blank_axis',color='wt1',/yst

vline,/h,2.35,thick=3,color='crimson'

cgPlots,432,2.24,psym='filleduptriangle',color='pur2',symsize=1.3
cgPlots,432,2.24,psym='openuptriangle',color='black',symsize=1.3
;; multiplot

;; ;; Panel 4
;; cgPlot,iter.nfit,iter.alpha,psym='filledcircle',xtit='N!d>min!n',$
;;        charsize=1.0,yr=[1.5,4.5],/yst,symsize=0.7
;; cgPlots,95,1.81,psym='filleduptriangle',color='bisque',symsize=1.3


myps,/done,/mp



END
